function h = generateColorHist(img, nbins, mode)
% GENERATECOLORHIST
% Input :
%   IMG   : Input color image.
%   NBINS : Number of bins in each color component.
%           Default is 8.
%   MODE  : case insensitive.
%      'NH'(default) for normalized histogram.
%      'H'           for histogram.
% Output :
%   H : Histogram, a nbins^3 matrix, in color component order.
%   For example, if input image is in RGB, K = 256/NBINS, histogram is:
%     R ~ K/2, G ~ K/2, B ~ K/2     -> H[(0*NBINS^2)+(0*NBINS)+1]
%     R ~ K/2, G ~ K/2, B ~ 3K/2    -> H[(0*NBINS^2)+(0*NBINS)+2]
%                    ...
%     R ~ K/2, G ~ 3K/2, B ~ K/2    -> H[(0*NBINS^2)+NBINS+1]
%     R ~ K/2, G ~ 3K/2, B ~ 3K/2   -> H[(0*NBINS^2)+NBINS+2]
%                    ...
%     R ~ 3K/2, G ~ K/2, B ~ K/2    -> H[NBINS^2+(0*NBINS)+1]
%     R ~ 3K/2, G ~ K/2, B ~ 3K/2   -> H[NBINS^2+(0*NBINS)+2]
%
%  P.S.: For greyscale image, gray value C is counted as (C, C, C).
%
img = double(img);
[P, Q, dim] = size(img);
if (dim ~= 1) && (dim ~= 3)
    error('Wrong number of color component, Use RGB_888(3) or Greyscale(1).');
end
if (max(img(:)) > 255) || (min(img(:)) < 0)
    error('Color component value must be in range 0..255.');
end
if nargin < 2
    nbins = 8;
end
if nargin < 3
    mode = 'nh';
end
k = 256 / nbins;
if (dim == 1)
    q = floor(img(:,:) / k)*(nbins*nbins+nbins+1);
else
    r = floor(img(:,:,1) / k)*nbins*nbins;
    g = floor(img(:,:,2) / k)*nbins;
    b = floor(img(:,:,3) / k);
    q = r+g+b;
end
h = histc(q(:), 0:nbins^3-1)';

if strcmpi(mode, 'nh')
    total = sum(h);
    if total > 0
        h = h / total;
    end
elseif strcmpi(mode, 'h')
    % do nothing
else
    error('Unknown mode.');
end
end
